Loading packages, creating seurat object

Set paths

# Set parent dir for visium experiments containing folders of spaceranger output with names
visium.dir <- "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RRRM2_RR10_test/RRRM2_ISST/Hippocampus"

# Path for where the plots will be exported to
export_path = "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RRRM2_RR10_test/RRRM2_ISST"

Create infoTable and create Seurat object

samples <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'filtered_feature_bc_matrix.h5')
spotfiles <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'tissue_positions.csv')
imgs <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'tissue_hires_image.png')
json <- list.files(visium.dir, recursive = TRUE, full.names = TRUE, pattern = 'scalefactors_json.json')
section.name <- samples
section.name <- gsub(paste0(visium.dir, "/"),"", gsub("/filtered_feature_bc_matrix.h5","",section.name))
infoTable <- data.frame(section.name, samples, spotfiles, imgs, json, stringsAsFactors = FALSE)

Data Quality Control

Number of Genes per Spot

VlnPlot(se_brain, features = "nFeature_RNA", ncol = 1, group.by = "section.name", pt.size = 0)

ST.FeaturePlot(se_brain, features = "nFeature_RNA", dark.theme = F, cols = c("lightyellow", "red", "dark red"), pt.size = 1.2, ncol = 4)

Number of UMIs per Spot

VlnPlot(se_brain, features = "nCount_RNA", ncol = 1, group.by = "section.name", pt.size = 0)

ST.FeaturePlot(se_brain, features = "nCount_RNA", dark.theme = F, cols = c("lightyellow", "red", "dark red"), pt.size = 1.2, ncol = 4)

VlnPlot - UMI/Genes/mtGenes/RiboGenes

se_brain$percent.mt <- PercentageFeatureSet(se_brain, pattern = "^mt")
se_brain$percent.ribo <- PercentageFeatureSet(se_brain, pattern = "^Rpl|^Rps")
VlnPlot(se_brain, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.ribo"), ncol = 2, pt.size = 0)

Run PCA here on the non-filtered data and then plot it

Ploting of the non-filtered data

p1 <- DimPlot(se_brain_merged_notfilt, group.by = "section.name")+ ggtitle("Plot of the non-filtered data")
p2 <- DimPlot(se_brain_merged_notfilt)

p1|p2

coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A","#88CCEE","#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499", "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")

DimPlot(se_brain_merged_notfilt, cols = coldef)

Filtering step

enids <- read.table(file = "/Users/solenefrapard/Documents/RRRM2_RR10_NASA_Project2/DataProcessing/RR3_test/rr3_ref/mm10_genes.tsv", header = T, stringsAsFactors = T)
enids <- data.frame(apply(enids, 2, as.character), stringsAsFactors = F)
rownames(enids) <- enids$gene_id
se_brain <- SubsetSTData(se_brain, expression = nFeature_RNA > 200)

mitochondrial, ribosomal genes and Malat1 filtered out

export Staffli object

staffli <- se_filtered@tools$Staffli
se_brain_merged <- merge(se_section_list[[1]], y = c(se_section_list[2:length(se_section_list)]), merge.data = T)
VariableFeatures(se_brain_merged) <- features

Plot data before and after harmony

PCA + Harmony

PCA

se_brain_merged <- RunPCA(se_brain_merged)
## PC_ 1 
## Positive:  Nrgn, Camk2n1, Slc17a7, Olfm1, Nptxr, Snap25, Cck, Chn1, Ppp3ca, Hpca 
##     Ctxn1, Calm2, Atp1a1, Rtn1, Snca, Calm1, Lamp5, Camk2a, Enc1, Ppp3r1 
##     Gpm6a, Ncdn, Mef2c, 1110008P14Rik, Aldoa, Hpcal4, Ptk2b, Cpne6, Stx1a, Slc30a3 
## Negative:  Plp1, Mbp, Mobp, Ptgds, Apod, Cnp, Trf, Mal, Fth1, Mag 
##     Cryab, Qdpr, Plekhb1, Scd2, Cldn11, Car2, Mog, Sept4, Csrp1, Gatm 
##     Dbi, Apoe, Sparc, Ndrg1, Ppp1r14a, Dbndd2, Ermn, Bcas1, Pllp, Mt1 
## PC_ 2 
## Positive:  Pcp4, Prkcd, Ramp3, Tcf7l2, Tnnt1, Rora, Adarb1, Ntng1, Plekhg1, Uchl1 
##     Cplx1, Pdp1, Ptpn4, Ndrg4, Slc17a6, Ccdc136, Plcb4, Zic1, Ptpn3, Amotl1 
##     Shox2, Rgs16, Nsmf, Pcp4l1, Patj, Kcnc2, Synpo2, Atp2a2, Nrip3, Gabra4 
## Negative:  Nrgn, Fth1, Ttr, Plp1, Cnp, Mal, Tmsb4x, Mobp, Enpp2, Apod 
##     Ddn, Clu, Ctxn1, Cryab, Trf, Gfap, Mbp, Aplp1, Cldn11, Hpcal4 
##     Cst3, Hpca, Camk2a, Penk, Mag, Ptgds, Camk2n1, Cplx2, Psd, Nptxr 
## PC_ 3 
## Positive:  Nrgn, Plp1, Mbp, Fth1, Snap25, Cck, Mobp, Slc17a7, Camk2n1, Pcp4 
##     Prkcd, Mal, Scn1b, Cnp, Mag, Cryab, Cplx1, Slc24a2, Trf, Hpca 
##     Chn1, Sept4, Rasgrp1, Ppp3ca, Adcy1, Qdpr, Car2, Dbndd2, 1110008P14Rik, Mog 
## Negative:  Nnat, Resp18, Hap1, Nap1l5, Ndn, Zcchc12, Calb2, Gaa, Gpx3, Sparc 
##     Tmem130, Ahi1, Baiap3, Ly6h, Wdr6, Gap43, Cartpt, AW551984, Scg2, Nrsn2 
##     Agt, Fxyd6, Gprasp2, Pnck, Pgrmc1, Vat1l, Tac1, Ngb, Bex1, Dlk1 
## PC_ 4 
## Positive:  Plp1, Mbp, Fth1, Resp18, Mobp, Nap1l5, Cnp, Stmn1, Mag, Mal 
##     Qdpr, Ndn, Stmn3, Sept4, Tubb4a, Cryab, Tmsb10, Gap43, Uchl1, Mog 
##     Tmem130, Trf, Ubb, Zwint, Dbndd2, Bex2, Zcchc12, Hap1, Gatm, Gapdh 
## Negative:  Apoe, Cst3, Clu, Slc1a2, Hbb-bs, Camk2a, Mt1, Atp1a2, Igfbp2, Glul 
##     Igf2, Hba-a1, Hba-a2, Ttr, Ddn, Mt2, Pltp, Ptgds, Lars2, Id3 
##     Mgp, Vtn, Rbp1, Mt3, Fxyd1, Hbb-bt, Vim, Gfap, Dbi, Ifitm3 
## PC_ 5 
## Positive:  Penk, Pcp4, Ncdn, Ppp1r1b, Gpr88, Pde1b, Hpca, Rgs9, Wipf3, Wfs1 
##     Gng7, Tac1, Ddn, Pde10a, Fkbp1a, Prkcd, Adora2a, Adcy5, Necab2, Ppp3ca 
##     Ptk2b, Prkcg, Camkv, Tmsb4x, Cpne6, 2010300C02Rik, Lrrc10b, Cnih2, Crym, Atp2b1 
## Negative:  Snap25, Pvalb, Camk2n1, Lamp5, 1110008P14Rik, Nrgn, Atp1a1, Stmn1, Kcnab3, Vsnl1 
##     Cabp1, Mef2c, Scn1b, Osbpl1a, Stx1a, Ttr, Serpini1, Cplx1, Nap1l5, Syt2 
##     Vamp1, Egr1, Igfbp6, Olfm2, Hbb-bs, Gabra1, Gnas, Sncb, Arc, Nrsn1
se_brain_merged <- RunUMAP(se_brain_merged, dims = 1:35, reduction = "pca")
## 09:49:38 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:49:38 Read 41301 rows and found 35 numeric columns
## 09:49:38 Using Annoy for neighbor search, n_neighbors = 30
## 09:49:38 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:49:41 Writing NN index file to temp file /var/folders/t8/1flxfkx50cb3ycdyf5y0gb0h0000gn/T//RtmpDtghuS/file12ee3ad9ebd8
## 09:49:41 Searching Annoy index using 1 thread, search_k = 3000
## 09:49:50 Annoy recall = 100%
## 09:49:50 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:49:52 Initializing from normalized Laplacian + noise (using irlba)
## 09:49:53 Commencing optimization for 200 epochs, with 1892064 positive edges
## 09:50:18 Optimization finished
p3 <- DimPlot(se_brain_merged, group.by = "section.name")+ ggtitle("Plot of the filtered data before integration with Harmony")
p4 <- DimPlot(se_brain_merged)

p3|p4

coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A","#88CCEE","#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499", "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")

DimPlot(se_brain_merged, cols = coldef)

Harmony

Clustering

p7 <- DimPlot(se_merged_harmony, group.by = "section.name")+ ggtitle("Plot of the filtered data after integration with Harmony")
p8 <- DimPlot(se_merged_harmony)

p7|p8

coldef <- c("#004A55","#0000FF", "#0031D4", "#00CBD4", "#006A21", "#799400","#FFD66A", "#D4C200", "#D6E279","#D4D6FF", "#D3AAFF", "#7500BF","#AA0084", "#FF9A94", "#FF8861","#D43938", "#946E55", "#7F1100", "#9B9A9A","#88CCEE","#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499", "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")

DimPlot(se_merged_harmony, cols = coldef)

se_merged_harmony@tools$Staffli <- staffli
ST.FeaturePlot(se_merged_harmony, features = "seurat_clusters", dark.theme = F, pt.size = 2.3, show.sb = F, ncol = 4, cols = coldef)

#top5_2 <- markers_RR10_2 %>%
top5_2 <- markers_RR10 %>%
  group_by(cluster) %>%
  filter(p_val_adj < 0.01) %>%
  filter(row_number() %in% 1:5)

d2 <- DotPlot(se_merged_harmony, features = unique(top5_2$gene) %>% rev()) + 
  coord_flip() + 
  scale_colour_gradientn(colours = RColorBrewer::brewer.pal(n = 11, name = "RdBu") %>% rev()) + labs(y = "cluster")
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
d2